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ABSTRACT 

The consistency of the Starburst model for AGN is tested using the optical variability 
observed in large data bases of QSOs. Theoretical predictions for the variability- 
luminosity relationship and structure function are presented and compared with ob- 
servations. If QSOs follow a variability-wavelength relation as that observed in nearby 
AGN, the model proves successful in reproducing the main characteristics of optical 
variability. The wavelength dependence (1) flattens the, otherwise, monochromatic 
Poissonian variability-luminosity relationship; and (2) decreases the asymptotic value 
of the structure function, which reveals that the elementary pulse driving the variations 
would have a characteristic time scale of 85-280 days. The upper limit is consistent 
with the time scale found in nearby Seyfert galaxies. Shorter values of this time scale 
are expected if the metallicity of high redshift objects is high, as recent observations 
indicate. If distant QSOs do not follow the variability-wavelength dependence ob- 
served in Seyfert nuclei and nearby QSOs, the characteristic pulse of variation needs 
to be much faster in order to reproduce the variability-luminosity relationship but, 
then, the single-parametric model explored in this work predicts a more rapidly rising 
structure function than that inferred from the data. 
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INTRODUCTION 



Variability is one of the most conspicuous properties of ac- 
tive galactic nuclei (AGN) and a potentially powerful con- 
straint on models for these objects. In the context of the 
Starburst model, the variability observed in Seyfert galax- 
ies and QSOs is produced by supernova explosions (SNe) 
which generate rapidly evolving compact supernova rem- 
nants (cSNRs) due to the interaction of their ejecta with 
the high density circumstellar environment created by their 
progenitor stars. During this phase, when the stellar cluster 
is 10-60 Myr old, the bolometric luminosity is dominated 
by stars, while the basic broad line region (BLR) proper- 
ties can be ascribed to the evolution of cSNRs in a medium 
with densities n <; 10 7 cm -3 and metallicities of the order 
of Zq or higher (Terlevich et al. 1992). In fact, the energy 
and overall pattern of variability of the optical light curves 
of intensively monitored Seyfert nuclei, such as NGC 4151 
and NGC 5548, can be well modeled by a sequence of SN 
and cSNR events (Aretxaga & Terlevich 1993, 1994) and the 
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detailed response of the BLR to these variations can also be 
explained by the evolution of the physical properties of the 
structures created by the remnants (Terlevich et al. 1995), 
with basically just one free parameter, the density of the 
ambient medium. Interestingly, the values obtained for this 
parameter in <; 10 7 cm -3 ) are also those derived from ob- 
served cSNRs, the so called "Seyfert 1 impostors". These 
cSNRs present optical spectra that resemble very closely 
those of type 1 Seyfert nuclei and QSOs (Filippenko 1989, 
Stathakis & Sadler 1991, Turatto et al. 1993; see Terlevich 
1994 for a review of their properties) . 

In the last decade, several observational studies have 
reported trends in the variability characteristics of optically 
selected quasars with either luminosity or redshift (Pica 
& Smith 1983, Cristiani, Vio & Andreani 1990, Giallongo, 
Trevese & Vagnetti 1991, Smith et al. 1993, Hawkins 1993, 
Hook et al. 1994; Trevese et al. 1994, Paltani & Courvoisier 
1994, Cristiani et al. 1996, Di Clemente et al. 1996). Some 
of these studies (see for example Pica & Smith 1983) have 
shown that the observed anti-correlation between optical 
variability and QSO luminosity is somewhat flatter than a 
"1/vA" law, a result which seems to rule out simple Pois- 
sonian models. However, these studies disregard the fact 
that AGN variability is wavelength dependent. In Seyfert 
nuclei and low-redshift QSOs it is observed that the ampli- 
tude of the variations increases towards shorter wavelengths 
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(Edelson, Krolik & Pike 1990, Kinney et al. 1991, Paltani 
& Courvoisier 1994, Di Clemente et al. 1996). Extrapolat- 
ing this behaviour to high-redshift QSOs, one expects the 
variability measured at a fixed optical band to overestimate 
the monochromatic rest frame optical variability, simply be- 
cause the objects are observed at bluer emitted wavelengths. 
Wavelength effects must be removed before analyzing the 
variability dependence with luminosity, and certainly be- 
fore claiming that the Poissonian model is inconsistent with 
the data. Indeed, the parametric fits to the wavelength- 
luminosity-redshift relationship of QSOs performed by Cid 
Fernandes, Aretxaga & Terlevich (1996) show that QSO 
variability can be Poissonian, once wavelength effects are 
considered. 

In this work we present the Poissonian model for the 
optical variability of QSOs as derived from the Starburst 
hypothesis, and compare it to the observations. The model 
makes specific predictions that are easy to check, namely, 
that the variability should be the result of a random su- 
perposition of events with luminosity and shape typical of 
cSNRs, and with a rate given by the cluster luminosity. The 
optical variability properties of massive stellar clusters con- 
taining cSNRs are analytically predicted in Section 2. In 
section 3 we use the South Galactic Pole sample of nearly 
300 QSOs to test the consistency of the theoretical predic- 
tions. An analysis of wavelength-variability dependence is 
presented there. Monte Carlo simulations are used to re- 
produce the sampling and photometric uncertainties of the 
observations. In Section 4 we discuss the results of the mod- 
elling of the variability-luminosity relationship and struc- 
ture function. Section 5 summarizes our main conclusions. 



2. THE VARIABILITY PRODUCED BY 

YOUNG CLUSTERS CONTAINING CSNRS 

The variability generated by a Poissonian process is charac- 
terized by the following properties: 

(i) the rate of events; 

(ii) the energy of the events; 

(iii) the shape of each variation pulse, i.e. its time evolu- 
tion; 

(iv) the strength of a non-variable "background" compo- 
nent, if it exists. 

The Starburst model makes specific predictions for each of 
these ingredients. In fact, we shall see that the model reduces 
the description of the variability properties of QSOs to a 
single functional parameter. 

Stars and cSNRs contribute to the mean B band lumi- 
nosity of a coeval stellar cluster in an approximately con- 
stant proportion during the SN II phase (10-60 Myr). This 
is so because both the SN rate (^sn) and the optical lumi- 
nosity coming from stars (L B ) are linked to the number of 
massive stars present in the cluster at every moment. The 
ratio of these quantities is given by 



2 x 10" 



(1) 



almost independently of the Initial Mass Function (IMF) 
and age of the cluster (Aretxaga & Terlevich 1994). The 
number of low mass stars, which carries the bulk of the 
mass of the cluster, is irrelevant in this ratio. The mean 



total luminosity of the cluster (Z/ B lu ) is given by the sum of 
the luminosity coming from stars and the mean luminosity 
coming from SNe (L B Nc ). From expression (hi) , this is related 
to the SN rate by 



Lclu T * I T SNe 
B — J^B T J^B 



5 x 10 10 i^L ( £B + i) Le 



(2) 



where eb is the mean B-band energy released in each SN 
remnant, in units of 10 51 erg. An estimation of ea can be 
obtained from the observed time averaged equivalent width 
of H/3 (Aretxaga & Terlevich 1994), 



Whp ~ 320A 



CB 



(3) 



l + 0.17e B 

This expression is also independent of the age, mass or IMF 
of the cluster, but is weakly dependent on the adopted bolo- 
metric correction for cSNRs. Most of the optical-UV lumi- 
nosity of cSNRs is emitted by dense shells of relatively cold 
gas formed during the remnant evolution, which reprocess 
a large fraction of the high energy photons generated by 
the fast shocks (Terlevich et al. 1992, Franco et al. 1994). 
The calibration of the bolometric correction was carried out 
using the intensities and equivalent widths of the observed 
cSNRs SN 1987F and SN 1988Z, giving e B /e 5 i w 0.12 where 
£51 is the energy emitted in B-band in units of 10 erg. How- 
ever, any systematic effect in the bolometric correction, due 
to metallicity or density differences, for example, could af- 
fect the energy estimation derived from expression (^) . The 
mean equivalent width of H/3 for typical QSOs is found to 
be Wh^ « 100 A (Searle & Sargent 1968, Yee 1980, Shuder 
1981, Osterbrock 1991), and from equation (^) the energy of 
cSNRs is derived to be £b ~ 0.5 (esi ~ 4). Values of the ki- 
netic energy released in a SN explosion in excess of 10 51 erg 
have been measured in recent well followed up type II SNe 
(Branch et al. 1981, Woosley 1988). 

The basic element of variation in the Starburst model 
is generated by a SN explosion and the evolution of its as- 
sociated cSNR. A simple description of such an event can 
be obtained through semi-analytical solutions (Shull 1980, 
Wheeler, Mazurek & Sivaramakrishnan 1980, Terlevich et 
al. 1992) as a double peak light curve given by (Aretxaga & 
Terlevich 1994) 



'(*) 
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1st peak: SN 



2nd peak: cSNR 



— ( — ^ 
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if t = 

if t = 110 days 

if t = 0.3 i sg 
if t > i S g 



Linear interpolation is used from the zero luminosity lev- 
els in both peaks. The zero point in the time scale corre- 
sponds to the SN outburst and t sg denotes the time beyond 
which radiative cooling becomes important in the evolution 
of the cSNR. For solar metallicity, t sg is given by (Shull 1980, 
Wheeler, Mazurek & Sivaramakrishnan 1980, Terlevich et al. 
1992) 



is 



0.62 yr e 51 1/8 n 7 



-3/4 



(5) 



where is the circumstellar density in which the cSNR 
evolves, in units of 10 7 cm -3 . This is a very schematic 
representation of the real luminosity evolution of cSNRs. 
Although the hydrodynamical simulations of Tenorio-Tagle 
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Figure 1. a. B band light curves of massive stellar clusters undergoing SN explosions with rates 0.1—500 yr - The rates are labelled 
inside each panel. The time of evolution of the cSNRs used for this set of models is t sg = 280 days and the mean B-band energy is 
eg = 0.5. The dashed line in the panels show the luminosity level produced by the stars of the cluster (following equation (1)). 



and co-workers (1996) show that the overall evolution of the 
luminosity of a cSNR roughly follows a i -11 / 7 law, several 
flares and fluctuations associated with shell formation, cool- 
ing instabilities and shell-shell collision occur on time scales 
shorter than t sg , with energies of about 10 49 -10 50 erg (Ter- 
levich el al. 1995, Plewa 1995). This general behaviour is 
also observed in isolated periods of activity in nearby type 1 
Seyfert nuclei such as NGC 4151, NGC 5548 and NGC 1566 
(Cid Fernandes, Terlevich & Aretxaga 1996, Alloin et al. 
1986). These secondary flares are not included in our sim- 



plified expression of the light curve. A further approximation 
made in equation (^) is that the B-band bolometric correc- 
tion remains constant throughout the lifetime of the cSNR. 

From these expressions we can deduce that the four 
properties that characterize a Poissonian process can be de- 
scribed with just one functional parameter in the Starburst 
model: 

(i) The rate of events is linked to the total luminosity of 
the objects by equation (Q) , and also depends on the energy 
of the events. 
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Figure l.b. As Figure l.a., for t sg = 85 days. 



(ii) The energy of the events in B-band is eb ~ 0.5 in 
order to satisfy equation Q with the observed values of the 
equivalent width of H/3 in QSOs. 



(iii) The main source of non-variable background lumi- 
nosity comes from the stars in the cluster. The stellar lumi- 
nosity is directly linked to the rate of events by equation . 
That makes for approximately 60 per cent of the total lumi- 
nosity of the nucleus. In R-band, the host galaxies of distant 
(z m 2) QSOs contribute to the non-variable light with less 
than 18 per cent of the total luminosity of the QSOs (Aretx- 
aga, Boyle & Terlevich 1995, Hutchings 1995), probably less 



in B-band. In what follows, the host galaxy contribution will 
be disregarded. 

(iv) The shape of the events is given by the SN+cSNR 
light curve of equation (^), which depends on the en- 
ergy (eb ~ 0.5) and the characteristic time scale of the 
pulse (t sg ). This last parameter actually controls the shape 
of the light curve, and remains free. However, its value 
can be constrained to a narrow band. For Seyfert nu- 
clei, such as NGC 4151 and NGC 5548, the values of t sg 
found to reproduce well isolated peaks in their light curves 
are 260-280 days (Aretxaga & Terlevich 1993,1994) but, 
since high luminosity QSOs may have higher metallici- 
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Figure I.e. As Figure l.a., for t sg = 15 days. 



ties (Hamann & Ferland 1992,1993), the evolution of their 
cSNRs could be substantially faster, as the cooling rates 
would increase. To allow for this effect we have set a lower 
limit of 15 days for i sg . 

Figures la,b,c show how the light curves of objects with 
different luminosities (Mb = — 19 to — 29 mag) look like, 
depending on the adopted value of the free parameter i ag , set 
to be 280, 85 or 15 days in our grid of models. To construct 
these model light curves of AGN, SN+cSNR light curves of 
a given t sg were combined at random under a certain rate, 
using a random number generator to determine the time 
of explosion. The final luminosity curve includes, in a self 



consistent way, the stellar luminosity as defined by relation 
(^). The light curves are sampled every 15 days, totaling 
100 years. SNe in a given cluster may have slightly different 
intrinsic parameters. For instance, their circumstellar densi- 
ties and/or B-band energies do not need to be exactly iden- 
tical. To allow for this possibility, we have let the values of 
tag and £b vary in a Gaussian way around their mean value, 
such that factor of 2 variations occur within a 2a probability. 
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2.1. Variability versus luminosity 

The light curve statistics of one of these clusters can be 
analytically calculated taking into account the Poissonian 
nature of the SN events (Cid Fernandes 1995). The mean 

luminosity due to SNe and cSNRs alone is simply Lb" = 
vsnEb, where E B = 10 51 £b erg is the average B-band energy 
radiated by each SN+cSNR event. The relative standard 
deviation of the SNe+cSNRs component is given by 



V B 



a(L B N °) 



(z'sntb 



,1/2 



(6) 



with tb being the effective lifetime of a cSNR in B-band, 
defined by 



j{f[L B ™(t;x)] 2 dt}p(x)di 



(7) 



where x denotes a given combination of t sg and e B and p(x) 
is the probability density of these two parameters. For a 
probability distribution as that of the simulations described 
in the previous section and for mean values of t sg and £b 
in the range 15 ^ t sg <, 280 days and 0.1 ^ £b ^ 0.5, the 
effective lifetime is found to be proportional to i sg , tb ~ 
5i sg . This scaling is almost independent of e B , i sg or the 
probability distribution of these two parameters. 

Besides the variable component due to SN events, the 
light curve also contains a stellar contribution which dilutes 
the variability (see equation (|l|)). The relative rms variability 
of the cluster, then, becomes 



„.clu 



C-B 



0.8 + ffl 



V B 



Using the scaling law between tb and t s 
ten as 



(8) 



can be rewrit- 



= 1.28e B 



^(l + 1.28e B ; 

yr 



The standard deviation v B lu is, thus, proportional to the in- 
verse square root of the average luminosity, as for a purely 
Poissonian process in which all the light is produced by in- 
dividual pulses. Note that equation (p|) is simply a "1 / y~N" 
law if all the pulses are identical (in energy and time-scale), 
but if there is any luminosity evolution in the pulse prop- 
erties it no longer describes a simple Poissonian process for 
the whole QSO luminosity range, although it is still ran- 
dom in nature. Equation ^ represents the same law, but 
diluted by the constant stellar background, which accounts 
for 60 per cent of the mean B-band luminosity of the cluster 
for cb ~ 0.5, as assumed here. The fact that the scaling with 
luminosity is the same in both equations is a direct conse- 
quence of the constant proportionality between the stellar 
luminosity and the SN rate (equation 1). 

Variability studies of optical light curves are usually car- 
ried out in magnitudes instead of luminosities. A problem 
arises when comparing published variability results with the- 
oretical predictions in that the logarithmic nature of mag- 
nitudes prevents the derivation of analytical predictions, 
as the ones presented above. One can, nevertheless, ob- 
tain approximate expressions which apply to the limiting 
cases of large and small v B u . In the limit of small v B u , 
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10 10 L® 



(9) 



t ag = 280 days 
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t =15 days 




M B (mag) 

Figure 2. Theoretical variability— luminosity relationship from 
the analytical predictions of v^ u and cr(M^ u ) for different t sg 
values. The abscissa in the upper panel is the absolute magnitude 
corresponding to L c ^ u , not the mean absolute magnitude of the 
cluster, as in the lower panel. 



i.e. high luminosity and SN rate, a small change in lumi- 
nosity can be transformed into a change in magnitude by 
AM ~ —(2.5 log e)AL/L. The rms variability in magnitudes 
can then be expressed as a(M B u ) ~ 1.09w B lu . However, for 
low luminosity objects, which according to equation ([]) have 
large v B lu , the first order expansion fails badly. In fact, it can 
be shown that a(M B lu ) increases proportionally to ^sn 1 ^ 2 in 
the limit of ^sntb <C 1. This property is most easily derived 
substituting our expression of L B SN (i) by square pulses last- 
ing tb and containing the same energy. Assuming that the 
pulses do not overlap in time we obtain 



a(M B 



2.5 [^sntb(1 



2.5 (fsN"r B ) 1/2 log 



™r 2 iog(i§ff +,s 



' 1.286b ' 
V^SNTB, 



(10) 



Figure 2 illustrates the variability-luminosity relation- 
ship for both v B lu and a(M B u ) and for different combinations 
of t a g. The two segments in the a(M B ) curves correspond to 
z'sntb < 0.5 (low luminosity) and i^sntb > 2 (high luminos- 
ity) regimes. 

These theoretical predictions for the rms variability 
must be used with caution when compared with observa- 
tions. While all the above expressions have been derived for 
well sampled, infinitely long light curves, observed QSO light 
curves usually span only a few years and are scarcely sam- 
pled. To illustrate these deficiencies, we have calculated the 
statistics of the light curves of figure la (t sg =280 days) for 
different sampling patterns, and compared them with the 
above predictions. The results are shown in figure 3. The 
cases considered are: (A) 7 points spaced by 1 yr, and (B) 
50 points spaced by 1 month. The filled squares and trian- 
gles in figure 3 show the results obtained for 10 realizations 
of cases A and B respectively. The large open squares and 
triangles represent the mean value of the ten realizations. 
The poor sampling in case (A) introduces a large extrinsic 
scatter in the plot, particularly at low luminosities. Even in 
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Figure 3. Comparison of the analytical variability— luminosity 
relationship with the relationship obtained from discrete finite 
light curves of (A) 7 points spaced by 1 yr, represented by squares, 
and (B) 50 points spaced by 1 month, represented by triangles. 
The small filled symbols represent the r.m.s. for fO individual 
realizations of cases (A) and (B), and the big open ones represent 
the mean value of these 10 realizations. 



case (B), which has a very good sampling compared to most 
observational data sets, it is clear that one cannot determine 
ct(Mb) to better than ~ ±0.05 mag for a Mb = —24 mag 
QSO, even in the absence of photometric errors. This uncer- 
tainty is larger when faster t ag regimes are used. 

A more severe problem is the bias induced by the finite- 
ness of the light curve in the measured variability of low 
luminosity objects. Low luminosity objects have a high in- 
trinsic variability because it is produced by a small num- 
ber of events. However this is only true for theoretical light 
curves. For observed light curves of low luminosity objects 
a is underestimated due to the low chance of catching an 
event in the time covered by the finite light curve. Note 
that in figure 3 most of the estimated a of objects with 
Mb <; —21 mag lay below the theoretical prediction. This 
effect is more and more important for objects of lower and 
lower luminosities, since the rate of events, and therefore the 
probability of catching an outburst decreases accordingly. 

Therefore, the model predicts a variability-luminosity 
anti-correlation for luminous systems. The relationship is 
dependent on the free parameter t sg , with larger amplitude 
fluctuations for shorter t sg values. For a finite light curve 
there is an inversion of the variability-luminosity relation- 
ship in low luminosity systems (Mb w — 20 mag), which for 
larger values of t sg is shifted towards lower luminosities. The 
inversion is caused by the increasingly lower probability of 
sampling a SN event in the time span of the observations. 



2.2. Variability as a function of time-scale: the 
structure function 

In recent studies of QSOs much effort has been concentrated 
on establishing the growth of variability with time, using 



the structure function (Bonoli et al. 1979, Simonetti, Cordes 
& Heeschen 1985, Cristiani, Vio & Andreani 1990, Hook 
et al. 1994, Trevese et al. 1994, Cristiani et al. 1996). The 
structure function (SF) of a light curve running from t = 
to T bs is defined as 

SF(r) = / [Lt(t + r)-L^(t)] 2 dt , (11) 

Tobs J t=Q 

and measures the mean squared luminosity variations of 
points spaced by a lag r. The SF is related to the auto- 
correlation function (ACF), 

ACF(r) = / L^(t + r)L^(t)dt , (12) 

fobs J t=Q 

by SF(r) = 2[ACF(0) - ACF(r)]. As in the case of rms 
variability, it is easy to derive analytical predictions for the 
behavior of the SF and ACF. The statistically expected ACF 
of the variable component, I/ B No (t), can be shown to be given 
by 



ACF (t) = N bN ACF (r) + 



(13) 



where N SN = ^sNiobs is the expected number of SN events 
in the light curve and ACF 1sn (t) is the ACF of an indi- 
vidual SN light curve, averaged over the distribution of its 
parameters, 



ACF (r) = 



1 

7~obs 



7"obs 



L B m (t + r- f)L B SN (i; x)dtp(x)dx 



(14) 



Note that the correlation between different events average 
out to the constant term in equation (p^), since the explo- 
sions may occur at any time in the light curve, with no pre- 
ferred time delay between events. Equation (^) does not in- 
clude the stellar background component, which simply adds 
an extra constant term to the ACF. The constant terms are 
canceled out in the SF, 



SF(r) 



= 2Ar SN [ACF 1SN (0) 
= JVSNSF^M . 



ACF 1 



(15) 



The shape of both the ACF and SF of the total light 
curve are thus determined by the luminosity profile of a 
single event. Since L B SN (t) tends to zero as the SN ages, 
ACF 1SN (r) tends to zero for large t, whereas the SF ap- 
proaches its asymptotic value 



SF(r 



co) = 2iV SN ACF 1SN (0) = 2a 2 (Ll Nc ) 



(16) 



The right hand side expression is obtained comparing the 
definition of ACF 1SN (0) with equations (g) and (0). The 
asymptotic limit of the SF is, thus, simply twice the pre- 
dicted variance of the light curve. This limit provides a nat- 
ural normalization which is very useful when analyzing the 
combined SF of many objects. The shape of the SF can- 
not be derived analytically for the SN light curve given by 
equation (^), but it is easily computed numerically. Figure 4 
shows the resulting SF normalized to its asymptotic value, 
plotted against the lag r in units of t ag . Plotted in this way, 
the SF is almost identical for different choices of es and t sg . 
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Numerical solution 
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Figure 4. Square root of the theoretical structure function (solid 
line) . Note that the representation is such that it is scale invariant 
with the parameter t sg . The dotted line represents the parame- 
terization of the SF given by equation (17). 



Figure 5. Luminosity— rcdshift plane for the QSOs in the SGP 
sample. Each QSO is marked by a dot. The luminosity bands 
marked by dashed lines enclose QSOs considered to have approx- 
imately equal SN rates, labelled inside the zone. 



This is due to the self-similar nature of L B SN (t) in the cSNR 
phase. From figure 4 we see that the SF of a starburst pow- 
ered QSO is expected to reach a flat regime for lags greater 
than 2-3t sg . Though an analytical result for the SF cannot 
be obtained, an acceptable parameterization of the numeri- 
cal result is given by 



SF(t) « SF(oo) 



1 + f 

t-sec 



(17) 



which is represented as a dotted line in figure 4. As for 
the rms variability, the SF is usually computed in terms 
of magnitudes, not luminosities. In the limit of small am- 
plitude variations, the relationship between these two SFs 
is similar to that between a(M^ u ) and u B lu , SFm ~ 

(2.51oge) 2 SF i /Zf :2 



3. COMPARISON WITH OBSERVATIONS 

Much of our knowledge about QSO variability is based 
on studies of large samples of sources monitored on pho- 
tographic plates over a period of one or two decades. In 
these kinds of studies we can place the Rosemary Hill Ob- 
servatory sample (McGimsey et al. 1975, Scott et al. 1976, 
Pica & Smith 1983, Pica et al. 1988, Smith et al. 1993), 
the ESO/SERC field 287 sample (Hawkins 1983,1986,1993), 
the Selected Area 57 sample (Trevese et al. 1989,1994) the 
Selected Area 94 sample (Cristiani, Vio & Andreani 1990, 
Cristiani et al. 1996) and the South Galactic Pole sample 
(Hook et al. 1994). 

We have chosen the South Galactic Pole sample (SGP) 
to test our models and to investigate the effect that the 
sampling of the sources, the photographic errors and the 
final ensemble of light curves might have in the deriva- 
tion of variability relationships. SGP comprises nearly 300 



QSOs, covering a wide band in the redshift-luminosity 
plane (0.3 < z < 4.1, -23 > M B > -29 mag for 
Ho = 50 Km s _1 Mpc -1 , qo = 0.5), observed for 7 epochs 
in a time span of 16 yr with relatively low observational er- 
rors (typically 0.07 mag). The distribution of QSOs in the 
luminosity-redshift plane can be seen in figure 5. 



3.1. The wavelength dependence of the observed 
variability 

Previous to the comparison of theory and observations, we 
have to consider that our predictions for light curves are 
computed in the rest frame B band, while real QSOs are ob- 
served in B band but their measured light has been emitted 
at shorter wavelengths. The mean magnitude of each QSO 
was k-corrected as in the original paper that presented and 
analyzed the data (Hook et al. 1994) in order to give the 
magnitudes represented in figure 5. In addition to this cor- 
rection, we must consider that the variability properties of 
the light emitted by QSOs from B band to up to A824 A, 
for the highest redshift QSOs of the sample (z = 4.1), are 
different if they behave similarly to low-redshift AGN. 

The UV studies of Edelson, Krolik & Pike (1990) and 
Kinney et al. (1991) show that the variability of nearby 
Seyfert nuclei and QSOs generally increases towards shorter 
wavelengths. The variability measured in these samples is 
between 40 and 60 per cent larger at A1450 A than at 
A2885 A. Similarly, for NGC 5548 we find that during the 
first AGN Watch campaign (Clavel et al. 1991, Peterson et 
al. 1992) the nucleus exhibited an intrinsic, i.e. photometric 
error free, cr(A1813 A) w 2.3 ctb, once the host galaxy con- 
tribution was subtracted, while for the whole 1978/90 his- 
tory of the nucleus (Wamsteker et al. 1989, Aretxaga 1993) 
a(A1350 A) w 2.7ct b . In the case of NGC 4151 its 1978/1990 
optical (Snijders 1991) and UV (Clavel et al. 1987,1990) 
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nuclear light curves give <r(A1813 A) » 3.4 <tb- The ratio 
of r.m.s. at these wavelengths changes considering different 
parts of the light curves, but we can regard variability in 
UV to be, at least, a factor of 2 larger than in B band. A 
recent study of the UV variability of nearby QSOs [z < 1.3) 
proposes a linear variability-wavelength relationship in the 
rest frame range A1200-3200 A, such that o\jL\ changes by 
(6.2 ± 4.3) per cent for every 1000A (Paltani & Courvoisier 
1994). Extrapolating this result to the optical, we also ob- 
tain a factor 2-3 for the <r(A1350 A)/a B ratio, as in the case 
of NGC 4151 and NGC 5548. However, their parameteri- 
zation cannot be universal for the whole optical-UV range, 
since this would predict negative values of <tb for QSOs ob- 
served at rest-frame wavelengths A ^ 1200 A (z <; 2.5). Di 
Clemente et al. (1996) show that PG QSOs (z <> 1) do show 
a wavelength-variability dependence, variability being fac- 
tors of 2 to 3 larger at UV than at R band, depending on 
the time-scale considered (see their figure 4). 

While the fact that low-redshift QSOs and other AGN 
have larger variability amplitudes at shorter wavelengths 
has been demonstrated by several authors (Edelson et al. 
1990, Kinney et al. 1991, Paltani & Courvoisier 1994, Di 
Clemente et al. 1996), the functional form of that depen- 
dence is unknown. We shall explore two parameterizations of 
the variability-wavelength relationship using the result that 
nearby AGN vary about a factor 2 or 3 more at A1350 A 
than at A4200 A. Consider first a linear law, as suggested 
by Paltani & Courvoisier (1994), but parameterized in ra- 
tios of variability, such that it reproduces factors of 2 to 3 
changes in <r(A1350 A)/ob, 

^-°- 35 ^^Ulc^A- 4 - 2 ) +1 • < 18 > 

where a goes from 2 to 3 and, secondly, a power law that 
reproduces factors 2 to 3 changes in o-(A1350 A) /(7b, 

^fy = (42001) ' (19) 

where a goes from 0.6 to 1. Cid Fernandes, Aretxaga & 
Terlevich (1996) have recently analyzed the variability- 
luminosity-redshift relation for the QSOs of the SGP 
sample, concluding that the variability-wavelength anti- 
correlation is present in the data. Their parametric fits sug- 
gest a power law index a between 0.1 and 1.2. 

3.2. Simulations of the SGP sample 

In order to accomplish realistic comparisons between mod- 
els and observations we have run Monte Carlo simulations 
which reproduce the conditions of the observations in our 
grid of model light curves. Simulations are essential not only 
to reproduce the effects of the poor sampling of the data 
but, also, to investigate the effects that the ensemble of in- 
dividual light curves has in the finally measured variability 
properties of QSO samples. 

Following relation (Q) , a SN rate can be assigned to each 
QSO in the sample. Figure 5 shows the luminosity bands 
considered for the assignations of the model light curves. 
Typical QSOs with luminosities between —24 and —27 mag 
have SN rates between ~ 5 and 100 yr _1 . Since we do not 
have a priori information about the time evolution of the 
postulated cSNRs which may be producing the variability in 



these QSOs, we will consider the three t sg values of the mod- 
els shown in figures la,b,c (280, 85 and 15 days) separately. 
For each QSO in the sample a set of points is randomly se- 
lected in the corresponding model light curve, such that it 
reproduces the time intervals of the observations in the rest 
frame of the QSO, taking into account the time dilution fac- 
tor (1 + z)^ 1 . For each of these points a corresponding error 
is generated under the actual error distribution function of 
the QSO brightness band. The error distribution functions 
are described in Hook et al. (1994). The process is repeated 
a hundred times in order to provide a set of simulated sam- 
ples from which we can extract results with statistical sig- 
nificance for models with different t sg values. 

3.3. The variability— luminosity relationship 

The r.m.s. <t(Mb) of each simulated QSO light curve is mea- 
sured, its intrinsic photometric error is subtracted in quadra- 
ture and, afterwards, the median of the r.m.s. of QSOs 
within a luminosity band of 0.6 mag is calculated for indi- 
vidually simulated samples. Figure 6 shows the variability- 
luminosity relationship found in the simulations. The crosses 
represent the median value of the r.m.s. found for the 100 dif- 
ferent simulated samples, for each t sg value. The thick solid 
line marks the median value of the averaged a {Mb) variabil- 
ity index over the 100 simulated samples, and the dashed 
line shows the analytical approximation represented in fig- 
ure 2. Note that when the variability is high (a <; 0.1 mag), 
this analytical approximation fails to reproduce the median 
value of the variability of the simulations. 

The open circles in figure 6 represent the median of the 
r.m.s. ct(Mb) of the observed QSOs, calculated as with the 
simulations. The numbers enclosed in the open circles indi- 
cate the number of QSOs used in each luminosity bin. Open 
squares and triangles represent the derived B-band variabil- 
ity found after applying the linear variability-wavelength 
correction described by equation (|l^), where squares corre- 
spond to a factor 2 in the ratio cr(A1350 A)/ctb, and triangles 
to a factor 3. Similarly, filled symbols represent the variabil- 
ity found in the data after applying the power-law correc- 
tion given by equation (|is|), in order to reproduce factors 2 
(squares) and 3 (triangles) in the ratio cr(A1350 A)/<tb- 

The agreement between models and data must be 
judged by the intersection of circles, triangles and squares 
with the bands of crosses, which show the range of <jb found 
for the 100 simulated data sets. The model which best re- 
produces the SGP data (circles), without any variability- 
wavelength correction, is that of fast intrinsic variations 
t sg = 15 days. Objects fainter than —24 mag can also be 
fit by a model with moderately fast intrinsic variations, be- 
tween t sg = 85 days and i sg = 15 days. The model with 
Seyfert 1 like variations, t sg = 280 days, can definitely be 
ruled out at the luminosities of interest, except, maybe, for 
QSOs in the last bin (Mb « —28.3 mag), which, never- 
theless, reach the error subtraction uncertainty limit. How- 
ever, when a variability-wavelength correction is applied 
(squares and triangles), the data is best described in all lu- 
minosity bands by models with slow Seyfert 1 like variations 
(t S g = 280 days) or with moderate variations (i sg = 85 days), 
depending on the choice of variability-wavelength correc- 
tion. Data corrected with <r(A1350 A)/<tb ~ 2 laws (squares) 
are best described by moderate variations, while data cor- 
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Figure 6. Variability— luminosity relationship found in the Monte Carlo simulations of the SGP sample, compared with the measured 
relationship for the data. Crosses represent the values found for 100 simulated samples, and the thick solid line represents the median value 
of the r.m.s. measured in those 100 simulated samples. The dashed line represents the theoretical approximation for variability derived 
in magnitudes in the low variability regime Open circles represent the r.m.s. of the observed QSOs in SGP. Open squares and triangles 
represent the same relationship for the the SGP QSOs after applying the variability-wavelength correction given by equation (18), such 
that a factor 2 in the c(A1350 A)/ as ratio is represented by squares and a factor 3 by triangles. Filled squares and triangles represent 
the rms of SGP QSOs after correcting the wavelength dependence of variability with a power-law of indices —0.6 and —1.0 respectively 
(factors 2 and 3 for the same ratio), as given by equation (19). The numbers enclosed by the open circles indicate the number of 
objects considered in the corresponding bin of luminosity to derive the variability indices: 14,25,34,43,50,54,29,18,4 from Mg ss —23.5 
to —28.3 mag. 



rected with er( A1350 A)/<tb ~ 3 laws (triangles) are equally 
well described by Seyfert 1 like and moderate variations. The 
fast variation models (t sg = 15 days) give a poor description 
of the variability-wavelength corrected data, except for the 
last two luminosity bins (Mb —27.5 mag). 

3.4. Structure function 

We also measure in the simulated samples the ensemble 
structure function, defined as 

SF(Ai) =< (M Bi ,k - M Bl ,lf >median , (20) 

where k, I are all the combination of epochs of the i-th QSO 
that have a time interval [t*. — ti\ comprised by At, and the 
median is calculated over the magnitude differences of all 
QSOs in the bin of interest. This ensemble structure function 
considers all the individual light curves of QSOs as a whole, 
since the poor sampling prevents the derivation of structure 
functions for individual sources. 

The crosses in figure 7 represent the square root of the 
structure functions measured in the 100 simulated samples 
corresponding to each t sg value. The time bins for comput- 



ing median values are set to be a year for pairs of epochs 
with At > 1 yr, and half a year for pair of epochs with 
At < 1 yr, in order to map the rising of the structure func- 
tion. The thick dashes in figure 7 mark the median values 
of the 100 realizations. The open circles represent the en- 
semble structure function measured in the data, where the 
number enclosed indicates the number of magnitude differ- 
ences per bin. As above, squares and triangles represent the 
structure function of the data, once these are corrected for 
the variability-wavelength dependence with the parameter- 
izations given by equations ( |lg| ) (open symbols) and ( |l9| ) 
(filled symbols), to reproduce factors of 2 (squares) or 3 
(triangles) in the <r(A1350 A)/<jb ratio. 

From the comparison of the intersection of crosses, 
squares, triangles and circles in figure 7, we deduce that the 
way in which the variability of the QSO sample is attained 
is well described by models involving t sg = 85 to 280 days, 
if one allows for a variability-wavelength correction as that 
observed in nearby AGN. If no variability-wavelength cor- 
rection is applied, the data cannot be described by a sin- 
gle parameter model. Although the short variation models 
(t sg = 15 days) reproduce the final variability level of the 
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Figure 7. Square root of the structure function measured in the Monte Carlo simulations of the SGP sample, compared with the 
measured function of the observed sample. All the QSO light curves are treated as an ensemble light curve. Crosses represent the 
individual SFs obtained for 100 simulated samples and thick solid dashes represent the median values of all the simulated samples. Open 
circles represent the SF of the observed sample. Open squares and triangles represent the SF same after applying a linear variability- 
wavelength correction, such that a factor 2 in the <r(A1350 A)/ob ratio is described by squares and a factor 3 by triangles. Filled squares 
and triangles represent the data after correcting the wavelength dependence of variability with a power— law of indices —0.6 and —1.0 
respectively (factors 2 and 3 for the same ratio). The numbers enclosed by the open circles give the total number of differences involved 
in the derivation of the structure function for the corresponding bin of time: 1725,670,590,878,1170,1063,642 from At 0.5 to 0.55 yr. 



data with no variability-wavelength correction, the first two 
time bins of the SF show that the growth curve of this model 
is quicker than that exhibited by the data. The quality of 
the data does not allow to probe the models to a higher 
degree. The most important difference among models with 
different t S g values is the time scale in which the flattening 
point of the SF is to appear. This flat regime is attained in 
time scales that the present data is unable to map in greater 
detail. In order to have a better description of the rise of the 
SF, one would need observations with time intervals of the 
order of, at least, months. 

We must note that the ensemble SF defined by equa- 
tion ( |2o| ) combines QSO light curves with different variabil- 
ity amplitudes. Since QSOs of different amplitudes of vari- 
ability, on average, correspond to different luminosities (fig- 
ure 6), and those QSOs happen to be at different redshifts 
(figure 5), then, different time bins in the SF will be dom- 
inated by a different population of QSOs, i.e. those which 
happen to have pairs of observations with rest frame time 
intervals included in that particular time bin. If the mixture 
of QSOs changes from bin to bin, as is normally the case, 
the level of variability mapped for each bin in the ensemble 
SF will be dominated by differences in the degree of vari- 
ability of different mixtures, and those variability differences 
will mask the genuine shape of the SF of individual QSOs. 
We can see that the scale invariant SF of figure 4 is not re- 
produced in the simulated samples (figure 7), just because 
the asymptotic value of the SF for each simulated QSO is 
different. Therefore, when the SFs of different data sets are 
compared, the ensemble SFs can easily differ just because 
the dominant population of QSOs in each bin has a differ- 
ent degree of variability due to (a) a different distribution 



of the QSOs in the luminosity-redshift plane or (b) a dif- 
ferent observing sampling pattern. Differences in the SF of 
different data sets have already been reported by Cristiani 
and collaborators (1996). A normalized SF defined by 



SF n (At) =< 



Mb 



Mb 



°*{Mb) 



> median 



(21) 



where a, (Mb) is the r.m.s. of the light curve of the i-th QSO, 
would be a more powerful tool to establish differences in the 
intrinsic variability growth of different samples and popula- 
tions of QSOs, getting rid, at the same time, of variability- 
wavelength dependences. From a theoretical point of view, 
the normalized SF is a more natural way of computing the 
ensemble SF, since the SFs of individual QSOs would be 
normalized to their asymptotic values, equation (|l6|), be- 
fore combining them. However, extracting the errors to the 
normalized SF is plagued with biases which arise in the dis- 
tribution of normalized magnitude differences. This natural 
way of extracting information about the variability growth 
curve of QSOs will be investigated in a subsequent paper. 



4. DISCUSSION 

4.1. QSOs as scaled up Seyfert 1 nuclei 

We have shown that if high-redshift QSOs have a vari- 
ability dependence with wavelength similar to that of 
nearby Seyfert nuclei and low-redshift (z < 1) QSOs, with 
cr(A1350 A)/ <7B ~ 3 ratios, their observed variability prop- 
erties can be reproduced by a sequence of pulses with time- 
scales t sg ~ 280 days. This value is also that found to repro- 
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duce the isolated long-term variations observed in the light 
curves of intensively monitored type 1 Seyfert nuclei (Aretx- 
aga & Terlevich 1993,1994). Furthermore, the modelling we 
have performed in the previous sections adopts a linear pro- 
portionality between the number of individual pulses and 
the luminosity of the object (relation ^|). These results im- 
ply that any model that considers QSOs as scaled versions 
of Seyfert nuclei can successfully reproduce the variability 
relationships found for radio-quiet AGN over a range of lu- 
minosities spanning, at least, 8 mag. In other words, our 
results show that the characteristics of the light curve of 
a QSO of —27 mag can be explained superposing the light 
curve of a Seyfert nucleus of —21 mag, like NGC 5548, about 
250 times. The starburst model predicts such a behaviour 
through equation (Q), and gives a natural explanation for 
the linear proportionality between the number of events and 
the global luminosity of the source: the link between the SN 
rate and the B band luminosity originated in the stars of 
the cluster (relation 111) . 



4.2. The masses of the stellar clusters and 
possible metallicity effects 

In case the variability-wavelength relationship that ap- 
plies to distant QSOs is flatter than a factor 3 in the 
a(A1350 A)/cr B ratio, we have shown that the model re- 
quires shorter values for the intrinsic time of evolution of 
the postulated cSNRs in order to reproduce the variability- 
luminosity anti-correlation and the ensemble SF of QSOs. 
Values of er( A1350 A)/<tb ~ 2-3 would allow models with 
pulses of t sg w 85-280 days to reproduce the observed vari- 
ability relationships. Therefore, it is possible that the vari- 
ability of distant QSOs is produced by shorter intrinsic vari- 
ations than those observed in Seyfert 1 nuclei. In this sub- 
section we would like to interpret such a possibility in the 
framework of the Starburst model. 

The high SN rates required to explain high luminosity 
objects imply high masses for the stellar clusters. For typi- 
cal QSOs with luminosities between —24 and —27 mag, the 
masses of the postulated coeval stellar clusters vary from 
2 x 10 10 to 3 x 10 11 M , for a solar neighbourhood Initial 
Mass Function (see figure 3a of Aretxaga & Terlevich 1994) . 
The mass estimation changes with the shape of the IMF and 
the age of the cluster. These values of the masses derived 
from the SN rates are well inside the hypothesis of Terlevich 
& Boyle (1993) that the QSO phenomenon might correspond 
to the formation of the cores of nowadays normal ellipti- 
cal galaxies at z ~ 2. The most massive elliptical galaxies 
have masses of up to 10 13 Mq (Faber et al. 1989, Terlevich 
1992). If that is the case, distant high luminosity QSOs are 
expected to have high metallicities, as high luminosity ellip- 
ticals do (Faber 1972,1973, Mould 1978). In fact, high metal- 
licities of, perhaps, up to 10 Zq have been inferred from the 
emission line ratios of high luminosity QSOs (Hamann & 
Ferland 1993). The result of having a higher metallicity envi- 
ronment is translated into having shorter times of evolution 
for cSNRs. The reason is two fold, first one expects that 
stars would create a higher density circumstellar environ- 
ment (Abbott, Bohlin & Savage 1982, Kudritzky, Pauldrach 
& Puis 1987) which is able to reprocess the kinetic energy 
into radiation more rapidly and, second, radiative cooling 



of the gas is more efficient at high metallicities, which also 
accelerates the evolution of the cSNR. 

From the analysis of the SGP data set we cannot clearly 
isolate metallicity effects, i.e. variations of i sg with luminos- 
ity, although some variation might be reflected in the du- 
ality of the success of the models with t ss ~ 280 - 85 days 
for the variability-wavelength corrected data and, certainly, 
in the comparison of these values with those derived for the 
light curves of Seyfert nuclei: 260-280 days (Aretxaga & Ter- 
levich 1993,1994). If indeed there is any evolution in i sg , the 
variability process will not be a simple Poissonian one, as 
explained in section 2.1 



5. CONCLUSIONS 

We have presented analytical predictions for the variability 
expected in a young stellar cluster with cSNR. We find that: 

(i) Luminous clusters are characterized by large SN rates, 
such that the amplitude of their variations is low due to 
the combined effect of the superposition of events and the 
dilution caused by the underlying stellar cluster light 

(ii) The unit of variability (a SN+cSNR event) has an 
effective lifetime of a few years. The maximum degree of 
variability in a cluster is attained at the time when the vari- 
ations produced by old SNe die out, i.e. when the baseline to 
measure the variability is larger than a few years. Therefore, 
variability versus time scale of variability relationships are 
expected to be flat at large baselines. 

We have performed numerical simulations to compare 
the model predictions with the results obtained from the 
analysis of large data bases of QSOs. Our Monte Carlo sim- 
ulations were parameterized to reproduce the time cover- 
age and observational errors of one of the best QSO data 
bases available to date, the South Galactic Pole sample. 
The models we present have only one functional free pa- 
rameter: the evolutionary time scale of the cSNRs in the 
stellar cluster (tag)- Typical clusters with luminosities be- 
tween —24 and —27 mag have a rate of SN explosions of 
5-100 yr _1 . We have shown that the r.m.s. index and en- 
semble SF of the SGP sample can be reproduced with mod- 
els in which the intrinsic variations are described by t sg 
values about 85-280 days, once a wavelength-variability 
correction typical of nearby AGN is applied to the sam- 
ple of QSOs. The values of the i sg parameter obtained for 
type 1 Seyfert nuclei (Aretxaga & Terlevich 1993,1994) and 
observed cSNRs (Terlevich 1994) are also around 280 days. 
From these simulations we conclude that, if the wavelength- 
variability correction of QSOs is between the explored values 
of <r(A1350 A)/<tb ~ 2-3, corresponding to those observed 
in low-redshift QSOs: 

• The variability-luminosity anti-correlation measured in 
QSOs is consistent with that expected from a simple Poisso- 
nian process, and can be naturally explained by the multiple 
superposition of SN events. Given that previous work has 
shown that the light curves of Seyfert galaxies can be well 
modeled with a sequence of SN+cSNR events, we concluded 
that QSO variability can be regarded as a scaled up version 
of Seyfert variability. 

• The shape of the SF is well matched by that predicted 
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by SN+cSNR events. The flattening point in the structure 
functions of QSOs is the key signature of the characteristic 
evolutionary time scale of the postulated cSNRs. It is essen- 
tial to concentrate the observational effort into mapping the 
rising of the SF in monthly time-scales. Such mapping is a 
powerful discriminator among different parameters for the 
models. 

• The SN rates derived from QSO luminosities can repro- 
duce the observed variability inside the restricted range of 
i sg values considered ( 85 < t ss < 280 days ). The masses of 
the clusters derived from these rates, are in good agreement 
with those of the young cores of elliptical galaxies postu- 
lated by Terlevich & Boyle (1993) to model the luminosity 
function of QSOs. 

• Although we cannot identify a clear trend of t S g values 
with luminosity, the data is consistent with some evolution 
in the characteristic lifetime of the postulated cSNRs, as ex- 
pected from a higher metallicity environment at higher red- 
shift. This would imply a simple Poissonian process for each 
QSO, but with different pulse characteristics for QSOs of dif- 
ferent luminosities, so that the ensemble variability doen't 
correspond to that of a simple Poissonian process any longer. 

• The random superposition of individual variations in 
luminous systems creates smooth structures in the light 
curves with time scales of several years. These smooth varia- 
tions are actually observed in moderately well sampled light 
curves of QSOs (Pica et al. 1988, Hawkins 1993, Maoz et al. 
1994) 

Clearly, a deeper study of these effects is needed 
through, possibly, the combination of large data bases as 
SGP, with similar or lower observational errors. Larger sam- 
ples, a more complete covering of the luminosity-redshift 
plane and better sampling of the light curves would cer- 
tainly allow a more stringent test of the models presented 
here. 

More generally, we have shown that the wavelength de- 
pendence of variability must be taken into account when an- 
alyzing the variability properties of QSOs observed in a fixed 
band. However, the functional form of such a dependence 
is still empirically unknown and future work should con- 
centrate on establishing an accurate wavelength-variability 
relationship from the UV to the optical. The forms used in 
this paper were just exploratory, taking into account the fact 
that nearby Seyfert nuclei and low-redshift QSOs exhibit a 
(j(A1350 A)/<tb ~ 2-3. A detailed observational study on 
multi-wavelength variability of high-redshift objects should 
address the question of whether such an approximation is 
valid for high-redshift objects. 

A final caveat is that in all our analysis we have disre- 
garded the possible contribution to the emitted light by the 
underlying galaxy. At high redshifts this contribution may 
become important. Future work should explore this possible 
contribution. 
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